package fem2;

import inf.jlinalg.BLAM;
import inf.math.LinAlg;

import java.util.ArrayList;

import math2.FunctionBasisOnRnToR;
import math2.FunctionRnToR;
import math2.Geometry;
import math2.LinAlgUtilities;
import math2.LinearCombinationRnToR;
import math2.MathUtilities;
import fem2.enu.MeshPartType;
import fem2.enu.NeighbourCriteria;
import fem2.enu.SortingCriteria;

/**
 * Abstraction class for element of mesh (i.e node, edge, face).
 * 
 * @author hbui
 * 
 */
public abstract class MeshPart implements Comparable<MeshPart> {

	protected int d = 1;
	protected int[] dofNumbers;
	protected double[] uHatLocal;
	private FunctionRnToR[] uhs;
	protected FunctionBasisOnRnToR basis;
	private int id;
	private MeshPartType mpType;
	private SortingCriteria sortingCriteria = SortingCriteria.SORT_ON_ID;
	// sorting criteria to sort MeshPart collection

	private Element e = null; // associated element with this mesh part

	private Geometry geom = null;

	/**
	 * @return the associated id of this mesh part
	 */
	public int getId() {
		return id;
	}

	/**
	 * set the id for mesh part
	 * 
	 * @param id
	 */
	public void setId(int id) {
		this.id = id;
	}

	/**
	 * @return the sortingCriteria
	 */
	public SortingCriteria getSortingCriteria() {
		return sortingCriteria;
	}

	/**
	 * @param sortingCriteria
	 *            the sortingCriteria to set
	 */
	public void setSortingCriteria(SortingCriteria sortingCriteria) {
		this.sortingCriteria = sortingCriteria;
	}

	/**
	 * @return the mesh part type
	 */
	public MeshPartType getType() {
		return mpType;
	}

	/**
	 * set the mesh part type
	 * 
	 * @param mpType
	 */
	protected void setType(MeshPartType mpType) {
		this.mpType = mpType;
	}

	/**
	 * set the geometry for this mesh part
	 * 
	 * @param geom
	 */
	public void setGeometry(Geometry geom) {
		this.geom = geom;
	}

	/**
	 * get the geometry associated with this mesh part
	 * 
	 * @return
	 */
	protected Geometry getGeometry() {
		return this.geom;
	}

	/**
	 * set element for this mesh part
	 * 
	 * @param e
	 */
	protected void setElement(Element e) {
		this.e = e;
	}

	/**
	 * get the associated element
	 * 
	 * @return
	 */
	public Element getElement() {
		return this.e;
	}

	/**
	 * return number of dofs associated with mesh part
	 * 
	 * @return length of dofNumbers
	 */
	public int countDofs() {
		return this.countModes() * this.countFields();
	}

	/**
	 * 
	 * @return number of fields associated with the problem
	 */
	public int countFields() {
		return this.d;
	}

	/**
	 * Count the number of shape functions associated with the mesh part.
	 * 
	 * @return the number of shape function
	 */
	public int countModes() {
		// return nodes.length;
		int nModes = 0;
		Node[] nodes = getNodes();
		for (int i = 0; i < nodes.length; i++) {
			nModes += nodes[i].countModes();
		}
		return nModes;
	}

	/**
	 * get the space dimension associated with the mesh part.
	 * 
	 * @return 2 if 2D, 3 if 3D, 1 if 1D
	 */
	public int getDimension() {
		return getNodes()[0].getPosition().length;
	}

	/**
	 * Enumerate the dofs for the mesh part.
	 * 
	 */
	public void enumerateDofs() {
		int nModes = countModes();
		this.dofNumbers = new int[this.d * nModes];
		Node[] nodes = getNodes();
		int cnt = 0;
		for (int i = 0; i < nodes.length; i++) {
			int[] nodeDofs = nodes[i].getDofNumbers();
			System.arraycopy(nodeDofs, 0, dofNumbers, cnt, nodeDofs.length);
			cnt += nodeDofs.length;
		}
	}

	/**
	 * Returns the dof array associated with the mesh part.
	 * 
	 * @return dofNumbers
	 */
	public final int[] getDofNumbers() {
		return this.dofNumbers;
	}

	/**
	 * Returns the array of nodes associated with the MeshPart
	 * 
	 * @return array of nodes
	 */
	public abstract Node[] getNodes();

	/**
	 * Get the local solution vector.
	 * 
	 * @return uLocal
	 */
	public final double[] getUHatLocal() {
		return this.uHatLocal;
	}

	/**
	 * get position vector of the mesh part
	 * 
	 * @return
	 */
	public double[] getPosition() {
		int dim = getDimension();
		double[] Xe = new double[dim * countModes()];
		for (int i = 0; i < countModes(); i++) {
			double[] p = getNodes()[i].getPosition();
			for (int j = 0; j < dim; j++) {
				Xe[dim * i + j] = p[j];
			}
		}
		return Xe;
	}

	/**
	 * Sets the basis for this mesh part
	 * 
	 * @param b
	 */
	public void setBasis(FunctionBasisOnRnToR b) {
		this.basis = b;
	}

	/**
	 * Gets the polynomial basis of the mesh part
	 * 
	 * @return polynomial basis of the mesh part
	 */
	public FunctionBasisOnRnToR getBasis() {
		return this.basis;
	}

	/**
	 * set number of fields associated with the problem
	 * 
	 * @param d
	 */
	public void setNumberOfFields(int d) {
		this.d = d;
	}

	/**
	 * set local degree of freedoms at time t.
	 * 
	 * @param t
	 *            time
	 * @param u
	 *            global displacement vector
	 */
	public synchronized final void setUHatGlobal(double t, double[] u) {
		int n = countDofs();
		int[] dofs = getDofNumbers();
		Node[] nodes = getNodes();

		uHatLocal = new double[n];
		for (int i = 0; i < n; i++) {
			int gPos = dofs[i];

			if (gPos != -1) {
				uHatLocal[i] = u[gPos];
			}
		}
		/*
		 * set displacement values to a node in case it's prescribed (for
		 * homogeneous constraint)
		 */
		for (int i = 0; i < nodes.length; i++) {
			Node nd = nodes[i];
			Constraint cc = nd.getConstraint();

			for (int j = 0; j < this.d; j++) {
				if (cc != null && !cc.isFree(j)) {
					this.uHatLocal[i * this.d * nd.countModes() + j] = cc.getValue(j, t);
				}
			}
		}
		updateSolution();
	}

	/**
	 * this function is invoked when computation need to be restarted on the
	 * element. i.e nonlinear analysis
	 */
	public void updateSolution() {
		int d = countFields();
		if (this.uhs == null) {
			this.uhs = new FunctionRnToR[d];
		}
		int n = this.basis.getDimension();
		for (int i = 0; i < d; i++) {
			// if (this.uhs[i] == null) {
			// LinearCombinationRnToR lc = new
			// LinearCombinationRnToR(this.basis);
			// double[] uHatLocal = this.getUHatLocal();
			// double[] coef = new double[n];
			// for (int j = 0; j < n; j++) {
			// coef[j] = uHatLocal[d * j + i];
			// }
			// lc.setCoefficients(coef);
			// this.uhs[i] = lc;
			// } else {
			// for (int j = 0; j < n; j++) {
			// double tmp = uHatLocal[d * j + i];
			// ((LinearCombinationRnToR) uhs[i]).setCoefficient(j, tmp);
			// }
			// }

			LinearCombinationRnToR lc = new LinearCombinationRnToR(this.basis);
			double[] uHatLocal = this.getUHatLocal();
			for (int j = 0; j < n; j++) {
				lc.setCoefficient(j, uHatLocal[d * j + i]);
			}
			this.uhs[i] = lc;

		}
	}

	/**
	 * @param i
	 * @return the solution associated to according field active on the mesh
	 *         part
	 */
	public FunctionRnToR getSolution(int i) {
		// if (this.uhs == null) {
		// this.uhs = new FunctionRnToR[this.countFields()];
		// }
		// if (this.uhs[idx] == null) {
		// LinearCombinationRnToR lc = new LinearCombinationRnToR(
		// this.polyBasis);
		// double[] uHatLocal = this.getUHatLocal();
		// int n = this.polyBasis.getDimension();
		// double[] coef = new double[n];
		// int numFields = this.countFields();
		// for (int i = 0; i < n; i++) {
		// coef[i] = uHatLocal[numFields * i + idx];
		// }
		// lc.setCoefficients(coef);
		// this.uhs[idx] = lc;
		// }
		return this.uhs[i];
	}

	/**
	 * Returns value of shape function at reference point
	 * 
	 * *
	 * 
	 * <pre>
	 * |	N1	N2	.	.	|
	 * </pre>
	 * 
	 * @param xi
	 * @return
	 */
	public double[] valueAt(double... xi) {
		int n = getBasis().getDimension();
		double[] N = new double[n];
		for (int i = 0; i < n; i++) {
			N[i] = getBasis().getBasisFunction(i).valueAt(xi);
		}
		return N;
	}

	/**
	 * Evaluate physical coodinates at the specified point
	 * 
	 * @param xi
	 * @return
	 */
	public double[] XAt(double... xi) {
		return getGeometry().valueAt(this, xi);
	}

	/**
	 * Evaluate reference coodinates at the specified physical point
	 * 
	 * @param xi
	 * @return
	 */
	public double[] XiAt(double... x) {
		int dim = x.length;
		double[] xi = new double[dim];

		double[] f = XAt(xi);
		for (int i = 0; i < dim; i++) {
			f[i] -= x[i];
		}

		// Debugger.watch("f = ", f);

		double norm_f = LinAlg.norm2(dim, f);

		int cnt = 0;
		while (norm_f > MathUtilities.tol6 && cnt++ < 50) {
			double[] A = new double[dim * dim];
			double[][] tmp = jacXAt(xi);

			// Debugger.watch("jac = ", tmp);

			for (int i = 0; i < dim; i++) {
				for (int j = 0; j < tmp.length; j++) {
					A[i * dim + j] = tmp[i][j];
				}
			}
			LinAlgUtilities.backSolve(BLAM.TRANSPOSE, A, f);

			for (int i = 0; i < dim; i++) {
				xi[i] -= f[i];
			}

			f = XAt(xi);
			for (int i = 0; i < dim; i++) {
				f[i] -= x[i];
			}

			// Debugger.watch("f = ", f);
			norm_f = LinAlg.norm2(dim, f);
		}

		// Debugger.watch("number of iterations to compute xi = ", cnt);
		if (cnt >= 50) {
			Debugger.watch("x = ", x);
			for (int i = 0; i < getNodes().length; i++) {
				Debugger.watch(getNodes()[i].toString());
			}
			throw new Error("iteration does not converge at mesh part " + getId());
		}

		return xi;

	}

	/**
	 * Evaluates the Jacobian
	 * 
	 * <pre>
	 * | dx1/dxi1 dx2/dxi1	.	|
	 * | dx1/dxi2 dx2/dxi2	.	|
	 * |		.				.			.	|
	 * </pre>
	 * 
	 * at the specified point.
	 * 
	 * @param xi
	 *            natural coordinate
	 * @return the Jacobian matrix of size dim x dim
	 */
	public double[][] jacXAt(double... xi) {
		return getGeometry().jacobianAt(this, xi);
	}

	/**
	 * Computes Jacobian of solutions vector at reference point.
	 * 
	 * <pre>
	 * | duh1/dxi1 duh2/dxi1	.	.	|
	 * | duh1/dxi2 duh2/dxi2	.	.	|
	 * |	.	.	.	.	|
	 * </pre>
	 * 
	 * @param xi
	 *            natural coordinate
	 * @return jacobian of solutions at natural coordinate
	 */
	public double[][] jacUAt(double... xi) {
		int dim = getDimension();
		int d = countFields();
		double[][] dUdxi = new double[dim][d];

		for (int i = 0; i < d; i++) {
			double[] g = this.getSolution(i).gradientAt(xi);

			for (int j = 0; j < dim; j++) {
				dUdxi[j][i] = g[j];
			}
		}
		return dUdxi;
	}

	/**
	 * /** Computes gradients of shape functions w.r.t physical coordinates in
	 * an array
	 * 
	 * <pre>
	 * | dN1/dX1  dN1/dX2	.	|
	 * | dN2/dX1  dN2/dX2	.	|
	 * |    .        .    	.	|
	 * |    .        .    	.	|
	 * |    .        .    	.	|
	 * </pre>
	 * 
	 * @param xi
	 *            vector of natural coordinates
	 * @return matrix of gradients
	 */
	public abstract double[][] gGradientsAt(double... xi);

	/**
	 * Computes the Jacobian of the approximate solution uh w.r.t physical
	 * coordinates. The gradients of the components form the column vectors,
	 * i.e.
	 * 
	 * <pre>
	 * | duh1/dX1  duh2/dX1	.	.	|
	 * | duh1/dX2  duh2/dX2	.	.	|
	 * |	.	.	.	.	|
	 * </pre>
	 * 
	 * @param xi
	 *            natural coordinate of physical coordinate
	 * @return Jacobian of solutions w.r.t physical coordinate
	 */
	public double[][] gJacUAt(double... xi) {
		int d = countFields();
		int dim = getDimension();
		double[][] jq = jacXAt(xi);
		double[][] ju = this.jacUAt(xi);
		double[][] gju = new double[dim][d];

		LinAlgUtilities.invertNxN(jq);
		LinAlg.multiply(false, dim, d, dim, 1.0, jq, ju, gju);
		return gju;
	}

	/**
	 * check positive jacobian
	 * 
	 * @return true if positive, false if negative
	 */
	public boolean checkJacobian() {
		double[] origin = new double[getDimension()];
		double[][] jac = jacXAt(origin);
		double dJ = LinAlgUtilities.determinantNxN(jac);
		if (Math.abs(dJ) < MathUtilities.tol14) {
			throw new Error("zero jacobian detected");
		}
		if (dJ < 0) {
			return false;
		} else {
			return true;
		}
	}

	/**
	 * compute center of gravity without weight of the mesh part
	 * 
	 * @return center of gravity (without weight)
	 */
	public final double[] computeCenter() {
		int dim = getDimension();
		double[] c = new double[dim];
		Node[] n = getNodes();
		for (int j = 0; j < dim; j++) {
			for (int i = 0; i < n.length; i++) {
				c[j] += n[i].getPosition()[j];
			}
			c[j] /= n.length;
		}
		return c;
	}

	/**
	 * 
	 * @return area if mesh part is face in 2d<br>
	 *         volume for block in 3d<br>
	 *         length for edge in 1d
	 */
	public final double computeVolume() {
		IntegrationManager integrator = IntegrationManager.getInstance();
		double[] xi = integrator.getIntegrationPoint(getType(), 1, 0);
		double w = integrator.getIntegrationWeight(getType(), 1, 0);
		double[][] jac = jacXAt(xi);
		return w * LinAlgUtilities.determinantNxN(jac);
	}

	/**
	 * check if a point lie inside or in the boundary of mesh part
	 * 
	 * @param p
	 * @return
	 */
	public boolean checkInside(double... p) {
		/*
		 * TODO
		 */
		throw new NotImplementedError();
	}

	/**
	 * @return all mesh part which is neighbours with this mesh part
	 */
	public ArrayList<MeshPart> getNeighbours(Mesh mesh, NeighbourCriteria criteria) {
		throw new Error("method not implemented");
	}

	@Override
	public String toString() {
		String str = "[" + Integer.toString(getId()) + ": ";
		Node[] nodes = getNodes();
		for (int i = 0; i < nodes.length; i++) {
			str += " " + nodes[i].getId();
		}
		return str;
	}

	@Override
	public int compareTo(MeshPart o) {
		if (this == o) {
			return 0;
		}
		if (!this.getClass().getName().equals(o.getClass().getName())) {
			throw new Error("can't compare ojects of different type of MeshPart");
		}
		if (getSortingCriteria() != o.getSortingCriteria()) {
			throw new Error("can't compare ojects of different sorting criteria");
		}
		if (getSortingCriteria() == SortingCriteria.SORT_ON_ID) {
			if (getId() > o.getId()) {
				return 1;
			} else if (getId() < o.getId()) {
				return -1;
			} else {
				throw new Error(
						"duplicated mesh part. Each mesh part of the same type must keep an unique id");
			}
		} else {
			throw new NotImplementedError();
		}
	}

	/**
	 * check if the mesh part contain the node
	 * 
	 * @param n
	 * @return the index of node in mesh part nodes if the mesh part contains
	 *         the node, else return -1
	 */
	public final int contains(Node n) {
		Node[] nodes = getNodes();
		for (int i = 0; i < nodes.length; i++) {
			if (nodes[i].equals(n)) {
				return i;
			}
		}
		return -1;
	}

	/**
	 * get neighbouring indices for each type of mesh part
	 * 
	 * @param type
	 * @param i
	 * @return
	 */
	public static int[] getNeighbouringIndex(MeshPartType type, int i) {

		switch (type) {
		case TRI3:
			switch (i) {
			case 0:
				return new int[] { 1, 2 };
			case 1:
				return new int[] { 0, 2 };
			case 2:
				return new int[] { 0, 1 };
			}
		case QUAD4:
			switch (i) {
			case 0:
				return new int[] { 1, 3 };
			case 1:
				return new int[] { 0, 2 };
			case 2:
				return new int[] { 1, 3 };
			case 3:
				return new int[] { 0, 2 };
			}
		case TRI6:
			switch (i) {
			case 0:
				return new int[] { 3, 5 };
			case 1:
				return new int[] { 3, 4 };
			case 2:
				return new int[] { 4, 5 };
			case 3:
				return new int[] { 0, 1 };
			case 4:
				return new int[] { 1, 2 };
			case 5:
				return new int[] { 0, 2 };
			}
		case QUAD8:
			switch (i) {
			case 0:
				return new int[] { 4, 7 };
			case 1:
				return new int[] { 4, 5 };
			case 2:
				return new int[] { 5, 6 };
			case 3:
				return new int[] { 6, 7 };
			case 4:
				return new int[] { 0, 1 };
			case 5:
				return new int[] { 1, 2 };
			case 6:
				return new int[] { 2, 3 };
			case 7:
				return new int[] { 0, 3 };
			}
			/*
			 * TODO: implement below
			 */
		case HEX20:
			break;
		case HEX8:
			break;
		case LINE2:
			break;
		case LINE3:
			break;
		case QUAD9:
			break;
		case TET10:
			break;
		case TET4:
			break;
		default:
			break;
		}
		throw new Error("wrong parameters");
	}

	/**
	 * get list of edges (define by index of node) for a mesh part type
	 * 
	 * @param type
	 * @return
	 */
	public static int[][] getListOfEdges(MeshPartType type) {

		switch (type) {
		case TRI3:
			return new int[][] { { 0, 1 }, { 1, 2 }, { 0, 2 } };
		case QUAD4:
			return new int[][] { { 0, 1 }, { 1, 2 }, { 2, 3 }, { 0, 3 } };
		case TRI6:
			return new int[][] { { 0, 3 }, { 3, 1 }, { 1, 4 }, { 4, 2 }, { 2, 5 }, { 5, 0 } };
		case QUAD8:
			return new int[][] { { 0, 4 }, { 4, 1 }, { 1, 5 }, { 5, 2 }, { 2, 6 }, { 6, 3 },
					{ 3, 7 }, { 7, 0 } };
			/*
			 * TODO: implement below
			 */
		case HEX20:
			break;
		case HEX8:
			break;
		case LINE2:
			break;
		case LINE3:
			break;
		case QUAD9:
			break;
		case TET10:
			break;
		case TET4:
			break;
		default:
			break;
		}
		throw new Error("wrong parameters");
	}

}
